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The Nice model of|Gomes at al.| (|2005[) suggests that the migration of the giant planets 



(-H caused a planetesimal clearing event which led to the Late Heavy Bombardment (LHB) 

at 880 Myr. Here we investigate the IR emission from the Kuiper belt during the 
I history of the Solar System as described by the Nice model. We describe a method 

O for easily converting the results of n-body planetesimal simulations into observational 

properties (assuming black-body grains and a single size distribution) and further 
modify this method to improve its realism (using realistic grain properties and a 
three-phase size distribution). We compare our results with observed debris discs and 
evaluate the plausibility of detecting an LHB-like process in extrasolar systems. Recent 
surveys have shown that 4% of stars exhibit 24 /im excess and 16% exhibit 70 fim 
excess. We show that the Solar System would have been amongst the brightest of 
these systems before the LHB at both 24 and 70 /im. We find a significant increase in 
24 /im emission during the LHB, which rapidly drops off and becomes undetectable 
within 30 Myr, whereas the 70 /im emission remains detectable until 360 Myr after the 
LHB. Comparison with the statistics of debris disc evolution shows that such depletion 
^5 events must be rare occurring around less than 12% of Sun-like stars and with this 

level of incidence we would expect approximately 1 of the 413 Sun-like, field stars so 
far detected to have a 24 fj,m excess to be currently going through an LHB. We also 
find that coUisional processes are important in the Solar System before the LHB and 
that parameters for weak Kuiper belt objects are inconsistent with the Nice model 
interpretation of the LHB. 



H Key words: solar system:general - Kuiper Belt - circumstellar matter - planetary 

systems 



1 INTRODUCTION 

Over the past couple of decades, an increasing number of 
stars have been found to be orbited by discs of planetesimals 
and dust known as debris discs. As more and more discs are 
discovered it becomes possible to start building up a picture 



of how these debris discs evolve over time (see Wyatt|2008 
for a review). Recent surveys (e.g. Hillenbrand et al.||2008 
Trilling et al. 2008] [Carpenter et al.||2009| ) have shown that 



the number of Sun-like stars that have been observed with 
24 /im emission (produced by hot dust) decreases with age, 
but the number of stars with 70 fim emission (produced by 
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cold dust) remains approximately constant with age. These 
observations generally agree with models suggesting that de- 
bris discs evolve in steady-state becoming collisionally de- 



pleted over time (Lohne et al. 20081, although there are a 



few exceptions that have much more hot dust than would 



be expected from these collisional arguments (Wyatt et al 



2007a I 



The Solar System has its own debris disc, with the 
majority of its mass concentrated in the asteroid belt and 
Kuiper belt. These belts correspond to the hot dust and cold 
dust seen around other stars, but our own disc is much less 



massive than these observed discs ( Moro-Martm et al. 2008 1 



Simulations of accretion in the Kuiper belt and the forma- 
tion of binary Kuiper belt objects (KBOs) suggests that the 
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original Kuiper belt must have been much more massive 



for the largest objects to form (e.g. Stern 1996a Chiang 



et al. 20071, which leads to the 'missing mass problem' of 
the Kuiper belt as this mass deficit cannot be explained by 
coUisional processes alone. 

One model that does explain the missing mass of the 
Kuiper belt - along with the orbits of the giant planets and 
various other details of the structure of the Solar System - 
is the Nice model (|Gomes et al.,|2005] Tsiganis et al.,^2005[ 



MorbideUi et al.|2005| [Levison et al.|2008a| ). The Nice model 
was designed to explain the current orbital elements of the 
outermost planets ( Tsiganis et al.|[2005 l. It is based on the 
idea that the gas giants formed much closer together. Due 
to interactions with the planetesimal disc, Saturn, Neptune 
and Uranus migrated outwards and Jupiter migrated slightly 
inwards. When Jupiter and Saturn crossed their 2:1 mean 
motion resonance (MMR) the system became temporarily 
destabilised, affecting the orbital elements of the gas gi- 
ants. As Neptune moved out into the Kuiper belt, it dy- 
namically excited the orbits of many of the KBOs, causing 
them to evolve on to cometary orbits and impact the ter- 
restrial planets and moons. As the planets' orbits evolved, 
secular resonance sweeping would have excited the orbits of 



many of the asteroids (Gomes 19971, thus also causing a 



bombardment of asteroids on the planets and moons of the 
inner Solar System. Hence, the Nice model also explains the 
Late Ifeavy Bombardment (LHB) of the Moon - a period 
of intense bombardment in which most of the craters on the 
Moon were formed, which occurred around 3.9 billion years 



ago (Tera et al. 19741, of which the latest impactors were 
most likely main belt asteroids ( jKring fc Cohen|2002{|Strom| 
et al.|2005| |. 

This period of intense bombardment and dynamical de- 
pletion of the Kuiper belt is likely to have had a significant 
effect on the observable properties of the debris disc of the 
Solar System. In this paper we investigate this effect by con- 
verting the distributions of planetesimal mass from the Nice 
model into distributions of emitting surface area to discover 
how the Solar System would have appeared to a distant ob- 
server during its history. Although the Solar System's debris 
disc has a number of components, for this paper we concen- 
trate on the changes to the Kuiper belt and how this would 
have affected the observable properties of the Solar System. 
In section[2]we describe the Nice model data and the simple 
analytical model applied to it, which uses the assumption of 
black-body grains and a single slope size distribution. In sec- 
tion [3] we look into relaxing the assumptions of black-body 
emission and a single slope size distribution to discover how 
a more realistic model changes our initial conclusions. Our 
final conclusions are given in section [4] 



2 MODELLING 

2.1 The Nice model data 

Our model is based on planetesimal data from one of the 
Nice model runs ( [Gomes et al.|[2b05[ ). In the Nice model, 
the simulation begins with effectively 10,000 particles, each 
representing 1.05 x 10~* M© of Kuiper belt objects. Any 
particles that reached a heliocentric distance of 1000 AU or 
evolved onto orbits with perihelion, g < 1 AU were removed 



from the simulation. The data covers a 1.2 Gyr time period 
starting at the time at which the gas disc dissipates. [Gomes| 
et al. ( 2005 1 ran 8 simulations with varying disc inner radius. 
The data used here is for the run with a disc inner edge of 
~15.5 AU, which places the LHB at 879 Myr, close to the 



■^700 Myr given by the analysis of Strom et al 



(20051. This 



is also the most realistic of the Gomes et al. ( 2005 1 runs 
since particles within ~15.3 AU have dynamical lifetimes 
shorter than the gas disc lifetime showing that they would 
have disappeared by the time the gas disc dissipates. This 
run starts with an initial disc mass of 35 and has 24 Mq 
at the time of the LHB. If the disc is less massive than this, 
Jupiter and Saturn do not cross their 2:1 MMR and there is 
no LHB. If the disc is more massive then the final separation 
of Jupiter and Saturn is much larger than it is today. 

2.2 Mass evolution 

Using the orbital elements of the Nice model particles we can 
then calculate the position of the particles at each time-step. 
For now it is the one-dimensional distribution that we are in- 
terested in and so the particles are separated into radius bins 
to allow us to determine the mass distribution in the system. 
Due to the small number of particles in the simulation, the 
evolution of the mass distribution appears stochastic in each 
1 Myr time-step. To smooth out the evolution and make it 
more realistic, two changes are made. Firstly, each particle 
is replaced by 10 particles, each of which is one tenth of 
the original mass, and these particles are spread uniformly 
around the orbit in mean anomaly to simulate the entire 
range of radii that particle would have passed through dur- 
ing this time period. By doing this we lose any resonant 
structures present in the data but increase the resolution 
of the model. Secondly, the mass distribution at each time- 
step is averaged over 5 time-steps (~5 Myr). Thus any values 
given in this paper at a specific time are actually averaged 
over 5 Myr. 

Figure [l] shows the surface density of the disc: just be- 
fore the LHB (at 873 Myr); during the LHB (at 881 Myr); 
and at the end of the simulation (at 1212 Myr). Hence- 
forth we refer to these epochs as pre-LHB, mid-LHB and 
post-LHB. Before the LHB occurs, most of the planetesi- 
mals are confined to a ring ~15 AU wide, centred at about 
26 AU. The surface density profile interior to this ring (over 
the range 8-19 AU) has a slope of i?2 9±o.4 outside the 
ring (over the range 38-106 AU) the slope is The 
mass surface density within the ring is 2 orders of magni- 
tude above this. At the onset of the LHB, a large number of 
the planetesimals are scattered from the belt both inwards 
and outwards spreading out the distribution of mass. Al- 
though planetesimals are scattered inwards before the LHB, 
at the onset of the LHB the rate at which planetesimals 
are being scattered inwards is much greater than the rate 
at which they are scattered back out. This results in a sur- 
face density profile with a leading slope of 7ji-5±o i between 
1 and 27 AU and a trailing slope of i?-4 8±o.i between 27 
and 106 AU at 881 Myr. By 1212 Myr the planetesimals are 
highly scattered giving a surface density profile with a lead- 
ing slope (over the range 9-38 AU) of /j3 ''='=0 '* and a trailing 
slop e (over the range 38- 106 AU) of 7?-2 '^±" \ 

[Levison et al.| ( [2006[ ) ran simulations of ecliptic comets 
to help them understand the orbit of the comet 2P/Encke. 
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Figure 1. Mass distribution before, during and after tiie LHB. 

Their data can be used to find the spatial distribution of 
comets which we might expect to be similar to our dis- 
tribution at the end of the Nice model. Between 9-38 AU 
the number density of comets is proportional to R^'^ and 
between 38-106 AU the number density is proportional to 
Our slope for the inner region is steeper than that 



R- 

found by Levison et al. ( 2006 1 and there are no particles 



within 9 AU whereas their comet model includes particles 
as far in as 0.1 AU. In part this difference is because the 
Nice model removes particles as soon as g < 1 AU therefore 
underestimating the number of comets in this region. The 
slopes for the outer region compare much better despite the 



fact that Levison et al. ( 2006 1 only include objects that have 



interacted with Neptune in their simulations suggesting that 
cometary dynamics is broadly similar, at least in terms of 
spatial distribution. 

The Nice model ends roughly 3 Gyr ago. For this work 
we would like to extrapolate the post-LHB evolution so that 
we can compare the predictions of the Nice model with 
current observations of the Solar System and compare our 
model with observations of extrasolar debris discs. To do this 
we need to find how the mass of the system will continue to 
evolve after the end of the LHB. 

At the end of the Nice model run there are 322 particles 
remaining. Many of these particles have left the confines of 
the Kuiper belt, becoming comets and scattered disc objects 
(SDOs), with only a few remaining as classical Kuiper belt 
objects (CKBOs). Here we define CKBOs as objects with 
g > 38 AU and 42 < a < 47 AU and assume that, since 
these objects are now on orbits that no longer bring them 
close to Neptune (or any of the other planets), they would 
be expected to remain trapped in the classical belt for the 



rest of the Solar System's lifetime (Levison et al. 2008a I 



In reality, the mass of the CKB will be decreased due to 
chaotic diffusion by resonances, small KBOs encountering 
large KBOs and collisions. However, this reduction in mass 



is a small fraction of the total mass ( Gomes et al. 2008 1 



From this definition we find that 3 out of the 322 particles 
represent CKBOs and that the mass of the classical Kuiper 
belt remains fixed at Mckb = 0.010 ± 0.006 M®, which 
compares favourably with recent observational estimates - 
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Figure 2. Total mass of Kuiper belt objects in the Nice model 
and extrapolated mass beyond the Nice model using equation [l] 
The scattered disc (dashed line) and classical belt (dot-dashed 
line) contributions to the total extrapolated mass are also shown. 



[Fuentes fc Holman|2008| - although this may be an under- 
estimate as the more detailed modelling of [Levison et al?] 
(2008a I, which includes dynamical processes not present in 
[Gomes et aL] ( |2005[ ), gives a final mass of 0.05-0.14 M®. 

As there are no major changes to the dynamical pro- 
cesses in the system after the LHB, we assume that the dy- 
namical losses affecting the rest of the particles (for which 
we use the term primordial scattered disc) remain the same 
and so the total mass will continue to decline. Therefore, we 
can extrapolate the mass evolution to the present day and 
into the future (figure[2|. From ~950 Myr onwards, the total 
mass (in Earth masses) as a function of time (in Myrs) can 
be fitted by the equation: 

+ AfcKB (1) 



(l + (t- 995)/280)2 



which puts the current total mass at 0.03 M® and the mass 
in the primordial scattered disc at 0.02 M®, which agrees 
well with the observations that set the current mass of the 
scattered disc to between 0.01-0.1 M® ( Gomes et al.|[2008 



and references therein). The evolution of the total mass of 
Kuiper belt objects in the system is shown in figure [2] which 
also shows the constant component of the CKB and the 
depleting component of the SD which combine to make up 
the total extrapolated mass. 



2.3 Converting mass to dust emission 

Each particle in the simulation represents a collection of 
Kuiper belt objects of many different sizes which are as- 
sumed to be affected by dynamical perturbations in the same 
way. By making some assumptions about the size distribu- 
tion of the objects (which will be considered in more detail 
in section 



3.1 1 we can work out the cross-sectional area of 



0.008-0.1 M® (Gladman et al.||2001 Bernstein et al. 2004 



dust corresponding to the mass in each radius bin and thus, 
the fiux emitted from these particles. 

First we assume that a collisional cascade is set-up 
quickly (from t=0) and so the planetesimals are in colli- 
sional equilibrium with a differential size distribution of the 
form n{D) oc D^~^'''^ where D is the diameter of the parti- 
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cles (in km) and Qd = 11/6 for an infinite coUisional cascade 
(Dohnanyi 19681. This size distribution is assumed to ap- 



ply from the largest objects of size Dc down to the smallest 
particles of size Dbi- Particles smaller than this are created 
in collisions between larger objects but are blown out of the 
system by radiation pressure on dynamical timescales and 
so contribute little to the size distribution. 

Taking = 11/6, Dc = 2000 km, L>bi = 2.2 ^m and 
the particle density, p = 1000 kg m~^, we find that the to- 
tal cross-sectional area in each radius bin (o-(_R) in AU^) is 
related to the total mass in each radius bin {AI{R) in M^) 
by: 



a{R)_ 
M{R) 



: 0.19 Au^ m; 



(2) 



As a first approximation we assume that the particles 
act like perfect black-bodies, except at submillimetre wave- 
lengths where an additional modification is applied. This is 
fine for large grains but does not work for small grains. Im- 
provements to this assumption will be considered in section 
[3] Using the black-body emission from the particles and the 
cross-sectional area worked out in equation |2] we can find 
the flux density (in Jy) measured by a distant observer: 



= ^2.35 X 10""cr(i?)B4A,rbb(i?))d^^X; 



Tbb(-R) = 278.3Ly*i? 



l/4p-l/2 



(3) 



(4) 



where is the Planck function (in units of Jy sr^^), which 
is dependent on the wavelength and temperature, d is the 
distance to the observer (in pc), L* is the luminosity of the 
star (in units of Lq), R is the distance between the parti- 
cle and the star (in AU) and X\ is a factor that accounts 
for the drop off in the emission spectrum beyond ~200 fim. 
Here we take Xa = 1 for A < 210 /iia and Xx = A/210 
for A > 210 nm to be consistent with submillimetre obser- 



vations of extrasolar debris discs (see [Wyatt et al. 2007b I 



The numerical coefficient in equation I3| arises because dif- 
ferent units are used for different parameters, an approach 
employed throughout the paper. 

By plotting the flux density against wavelength for the 
dust emission (see figure [3| we can see how the LHB causes 
the emission spectrum to change. Before the LHB, the emis- 
sion resembles a single temperature spectrum appropriate to 
the radius of the belt (i.e. 55 K at 26 AU). During the LHB, 
the spreading of mass from the belt (see figure [TJ means 
that the dust is emitting from a much broader range of tem- 
peratures and so the spectrum covers a broader range of 
wavelengths. Wavelengths as low as 7 nm now have a flux 
density > 10"'' Jy as opposed to just wavelengths 16 fj,m and 
longer in the pre-LHB phase. In particular, we see that the 
mid-IR flux is enhanced. After the LHB has occurred, the 
flux at all wavelengths rapidly decreases and the spectrum 
begins to resemble a single temperature black body once 
again but at a longer peak wavelength due to the increase 
in mean radius of the belt. 

The increase in mid-IR flux seen here during the LHB 
is only a lower limit, since objects with perihelion less than 
1 AU are removed from the simulation, which removes a lot 
of planetesimals that are scattered onto cometary orbits (as 
discussed in section [2^ . These would otherwise contribute 
to the emission via dust production and sublimation. We 
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Figure 3. SED before, during and after the LHB, as it would 
appear from 10 pc away. The thick line shows the solar photo- 
sphere. The thin lines show the excess emission at 873 Myr, 881 
Myr and 1212 Myr. The total emission spectrum observed would 
be the sum of the photosphere and excess. 



also note that the true mid-IR emission may be higher at 
all times as we have only considered the contribution of the 
Kuiper belt and have not included the asteroid belt. These 
effects will be investigated in more detail in future work. 

Flux density is dependent on distance from the observer 
to the star. To be able to compare different debris discs it 
is necessary to use a variable independent of distance. As 
stellar flux density (in Jy) is also oc d~^: 



F^, = 1.77B^{X,T^)L^T-'^d- 



(5) 



the excess ratio {F^/F^t,, also called the fractional excess) 
is one distance independent measure, but is dependent on 
wavelength. A variable independent of both distance and 
wavelength is the fractional luminosity, /, which measures 
the ratio of the excess luminosity (due to the dust) to the 
luminosity of the star: 



id 



Jfvav 



(6) 



where Ld is the luminosity of the dust and is the lumi- 
nosity of the star. 

Figure |4] shows how fractional luminosity varies with 
time for our model. Before the LHB event, the fractional 
luminosity shows that the planetesimal disc is in a quasi- 
steady state in that dynamical losses of KBOs are a rela- 
tively small fraction of the total mass. At the time of the 
LHB there is a slight but minimal increase in / due to the 
influx of comets. After the LHB, / rapidly decreases in a 
similar manner to Mtot (see figure [2|. Here we assume that 
the radial distribution of the mass remains the same from 
the end of the Nice model. In reality, since the mass of the 
CKBOs remains constant and it is only the SDOs that are 



being lost through dynamical processes (see section 2.2 1, the 



distribution of mass will again resemble a narrow belt sim- 
ilar to that present before the LHB but at a larger radius. 
As such this assumption probably overestimates the mid-IR 
fiux at late times because there would not be as much mass 
spread inwards as we are assuming. 

By extrapolating the fractional luminosity we find that 
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Figure 4. Fractional luminosity as a function of time. Inset shows 
small increase in fractional luminosity during the LHB. However, 
it should be noted that the peak is diminished due to the fact 
that we have averaged over 5 Myr (see section ['2.2[ |. The fractional 
luminosity at late times is an overestimate since we have ignored 
P-R drag and SW drag effects (see section [2.5.2| . 



this model gives the current value to be / = 2 x 10 ^ which 
is within the range lO"'^ — 10~^ suggested by the size distri- 
bution of KBOs ( |Backman et al.|1995| [Stern|1996b' l. 



2.4 Comparison with extrasolar debris discs 

Dust in debris discs emits most strongly in the infrared, as 
can be seen in figure [S] Since its launch in 2003, the Spitzer 
Space Telescope has been used to survey stars for infrared 
excesses. The Multiband Imaging Photometer for Spitzer 
(MIPS) makes observations of the stars at 24 fim and 70 /im 
which can then be compared to photospheric models to cal- 
culate if there is evidence for any excess emission which may 
be due to dust present in the system. Surveys using Spitzer 
are generally calibration limited which means that they can 
detect stars with a fractional excess (the ratio of flux from 
the dust to flux from the star at a given wavelength) above 
a given limit. At 24 /xm, the Formation and Evolution of 
Planetary Systems (FEPS) survey can make 3ct detections 
of excess down to a limit of -^24/^24* = 0.054 for the bright- 
est stars ( Carpenter et al^|2009 |. At 70 fim the limit is ap- 
proximately Fjo/Frot « 0.55, although observations of the 
more distant stars are sensitivity limited and so have not be 
observed down to this limit ( Wyatt|[2008 1. 



Carpenter et al. 



(20091 surveyed 314 stars and found 
that there is a decrease in 24 fim excess with age. They 
show that 15% of stars younger than 300 Myr have a 24 fim 
excess greater than 10.2% above the photosphere but this 
fraction goes down to 2.7% for older stars. By combining 
observations of both field stars and stars in open clusters 
and associations from the literature, Caspar et al. ( 2009 1 



also find that there is a decrease in the fraction of stars with 
24 fim excess with age, levelling off at a few percent for stars 



older than 1 Gyr. Trifling et al. ( 2008 1 found that 16% of 



F and G type stars have detectable debris discs at 70 fim 
from a sample of 225 stars. Although they show that the 
data could indicate a decrease in the fraction of stars with 
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Figure 5. Excess ratio versus time for 24^m (top) and 70^m 
(bottom). The solid line represents the emission from our model, 
assuming a single-slope size distribution with = 11/6 and 
black-body grains (c.f. figure [T2I 1. The asterisks are observed discs 
and the dashed line shows the approximate observational limit. 
The excess ratio at late times is an overestimate since we have 
ignored P-R drag effects (see section [2.5.2| . 



detectable excess with age, a constant excess fraction also 
adequately flts the data and there are currently too few ob- 
servations to distinguish between the two.^Hillcnbran d et al.| 
(2008 ) similarly find no apparent trend in the 70 fim excess 
fraction with age, however they do note that the maximum 
excess ratio at 70 fim does appear to decrease with age, 
which can be seen in figure [5] (bottom) . 

The evolution of the fractional excesses at 24 ^m and 
70 fim for our model are shown in figure [5] For comparison 
these plots also show 106 Sun-like stars (represented by as- 
terisks) for which excesses have been detected (Habing et al 
pool; Bcic hman etl il. 2006; Moor et al . 2006, Trilling et al 
2007; Hillenbrand ot al. 2008; Trilling eTaTpoMnCarpem 



|ter et al.||2009j ). 77 of these stars have observed excesses at 
70 fim and 53 of them have observed excesses at 24 fim. The 
dashed lines show the approximate limits of detectability. 

From the 24 ^m excess, we can see that the hot emis- 
sion from the model starts at -F24/-F24* = 0.5, which is high 
enough to make the Kuiper belt detectable at early times. 
This hot emission gradually decreases during the pre-LHB 
phase and then briefiy rises again during the LHB back to 
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its initial value (see inset of figure[5](top)) due to an increase 
in the mass closer to the Sun (see figure [T]). The Kuiper belt 
is also detectable at 70 fim at early times. The 70 ^m ex- 
cess remains in a quasi-steady state until the LIfB at which 
point it drops off sharply, but still remains detectable up to 
360 Myr after the LHB. 

It is possible that some of the observed systems may 
be going through a similar process and that some systems 
may be observed whilst in the middle of an LHB-like epoch, 



especially those systems described in Wyatt et al. (2007a I 
as having a mid-IR excess (from dust at a few AU) higher 
than expected for their age. However, it is unlikely to explain 
systems like H[D69830 which have been detected at 24 fim 
but not at 70 ^m since our results imply that, although the 
70 /im excess of a system does decrease from the time of 
the LHB onwards, the system should still be detectable at 
70 fj,Tn for a few hundred million years after the LHB. In 
other words, a system must have a significant cold disc at 
the same time as the hot disc to provide material for the hot 
disc. 

Figure |5] (bottom) shows that there are a large number 
of observed discs at late times, which clearly have not gone 



through an LHB. The fact that the [TriUing et al.| ( |2008 1 
results are consistent with the fraction of Sun-like stars with 
a detectable 70 fj,m excess remaining approximately constant 
with age at 16.4^2 g% shows that extrasolar LHB events 
must be rare. Although, the number of systems surveyed 
is still fairly low, we can still place an upper limit on the 
fraction of systems that may undergo an LHB event. To 
get this limit we start by assuming that if a star is born 
with a planetesimal belt that is detectable at 70 /im then 
it remains detectable unless a major planetesimal-clearing 
event, like a late heavy bombardment, takes place. For the 
fraction of stars born with a detectable planetesimal belt we 
assume that the value of leAtl'lfo from TrilUng et al.| ( [2008 1 
also applies at the youngest ages (<100 Myr). Although the 
Trilling sample is not focussed on young stars, not including 
any systems younger than lOOMyr, the results of Carpenter 
et al. (2008) are consistent with the distribution of fractional 
excesses remaining constant for all ages. Thus the lack of 
decline tells us that the fraction of stars starting with a 
70 /im excess that go through a planetesimal-clearing event 
is 0% with a Sa upper hmit of 3^2(2.9/16.4)2 = 75%. This 
gives a maximum of 12% of all Sun-like stars experiencing a 
Late Heavy Bombardment event. 

Since we might expect an LHB event to require the 
presence of giant planets it is encouraging to find that this 
fraction is not greater than the fraction of Sun-like stars 
inferred to have gas giants (planets with masses equal to or 
greater than Saturn) within 20 AU which Marcy et al. ( 2005 1 
estimate as 12%. If LHBs were common for stars with giant 
planets then the presence of debris would be expected to 
be anti-correlated with the presence of giant planets for old 
stars. Since this is not observed to be the case ( [Greaves et al.| 
2004p and several old stars are now known with both giant 



^ Note that a direct comparison of the fraction of stars detected 
in each of these surveys is not possible since stars in the different 
surveys were observed down to different levels of fractional excess, 
notably with higher detection thresholds for the young stars in 
the Carpenter survey which are typically at greater distance. 



planets and debris (see table 1 in Moro-Martfn et al.||2007 
where all the stars are at least 500 Myr), this is further 
evidence that the fraction of stars that undergo LHBs is 
<12% (assuming that Saturn mass planets within 20 AU 
are required for an LHB). 



Caspar et al. ( 2009 1 also estimate the fraction of Sun- 
like stars that go through an LHB event. They find a maxi- 
mum limit of 15-30% based on observations of 24 /im excess. 
This is much higher than our own limit as they have been 
less restrictive with their definition of an LHB event. They 
assume that any star that is observed to have a 24 /xm ex- 
cess must be going through an LHB, however we have shown 
that debris discs can be detectable at 24 in the pre-LHB 
phase (see figure [5] (top) ) . Furthermore, only a small num- 
ber of systems have a 24 ^m excess that is too high to be 
explained by collisional processing ( Wyatt et al.|2007a l. 

Observations in the submillimetre part of the spectrum 
offer a useful method for estimating the dust mass of debris 
discs and so are often used as another method of comparing 
debris discs. The dust mass, Mdust in M®, can be calculated 
Zuckerman|2001| ): 

26 X 10^° F,d\-^ B-^ (7) 

where is the mass absorption coefficient (in AU^ ^e^)- 
By combining this with equation [3] we find that: 



using (e 

Mdust ~ 



(8) 



There is a lot of uncertainty in the value of Hi, since it is 
dependant on the properties of the particles in the system. 
Here we adopt the value K850m™ — 45 A U'^ M^'^ to ease com- 



parison with values reported elsewhere (Najita & Williams 
2005 1. Figure[6]shows the submillimetre dust mass predicted 
by our model. The asterisks represent dust masses for Sun- 
like stars that have been observed to have an excess at 
850 /xm. Data for these 13 stars is taken from the literature 



(Wyatt et al.||2003 


Creaves et al.||2004 


2005 


Sheret et al. 


2004 Najita & Wi 


liams||2005 Wyatt et al.||2005| |Williams 



fc Andrews|2006[ [Greaves et aL|20'09 | 



Creaves et al.[ ([2004[) use COBE/FIRAS observations 



at 800 /xm to provide an upper limit to the dust mass of the 
Kuiper belt, which they find to be ~2 x 10~^ Mg. Our model 
implies that the current dust mass is ~3.1 xl0~^ M®. The 
discrepancy between our result and the observations may 
be due to Poynting-Robertson drag being neglected in our 
model, which can have the effect of reducing the amount of 
small dust in a debris disc as described in section 12.5.21 

2.5 Collisional lifetime 

Mass loss in the Nice model is entirely due to the dynamical 
evolution of the particles. For computational reasons, it was 
assumed that the particles only interacted with the planets 
and not with each other. In reality collisional processes are 
likely to have had some effect on mass loss in the system. 
In this section we investigate the effect of collisions and P-R 
drag on our simple model. 

In section |2.3| we assumed that the particles are in a 
collisional cascade. In a collisional cascade, objects of size D 
to D -I- dD are destroyed by collisions only to be replaced 
by fragments created by collisions of larger objects. From 



the equations of Wyatt et al. ( 1999 2007a I it can be shown 
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Figure 6. Sub-mm dust mass as a function of time. The solid 
line represents tlie emission from our model and the asterisks are 
observed discs around F, G and K stars. The dust mass at late 
times is an overestimate since we have ignored P-R drag effects 
(see section [2.5.2[ l. 

that the time between catastrophic coUisions (known as the 
coUision time-scale) for particles of size D in a belt at a 
mean distance Rm (in AU) from the star and with a width 
dr (in AU) is given by: 

^ (9) 



tc{D) 



Ctot = 



Rjfdr 
aiR) 



2[l + 1.25(e//)2]- 



fcc{D) 



M{R) 



(10) 



where tc is in years, M^, is the mass of the star in solar 
masses, atot is the total surface area and e and I are the 
mean of the eccentricities and mean of the inclinations (in 
radians) respectively. /cc(-D) is a factor determined by the 
fraction of the total cross-sectional area which is seen by 
a particle of size D as potentially causing a catastrophic 
collision and is given by: 

f..{D)= il + D/D'fa{D')dD' (11) 

where Dcc{D) is the smallest particle that can catastroph- 
ically destroy a particle of size D and a is the normalised 
cross-sectional area distribution in each diameter bin. Since 
we are assuming a single power law with gd > 5/3 this can 
be written as: 



fcciD) = 



3gd-5 f Dcc{D) 



-3qd 



3gd 



+ 



2D{DrAD) 



4-3<7d 



+ 



3gd -4 
D2(£),,(D)3-3<,d _ 



D, 



3-39di 



(12) 



3(jd — 3 

DcciD) = X^D for XcD > and Dcc(-D) = other- 
wise. The factor Xc can be calculated using the equation: 

= 1.3 X 10-'[Q^i?n.M-V(e,/)"Y'" 



(13) 



where is the dispersal threshold and f{e,I) is the ratio 
of relative velocity to Keplerian velocity. This is given by: 



fie, I) = \/l.25e2 -h72. 



(14) 



Here we use the value Qd = 200 J kg ^ since this value 
provides a good fit to the statistics of debris discs around 
A stars ( |Wyatt et aL]|2007b| ). This is an effective planetesi- 
mal strength that describes the dust mass loss rate from the 
planetesimal belt which is linked to medium sized planetesi- 
mals (e.g. Dc = 160 km in|Wyatt et al.|2007b). In reality 



varies with size. Thus we expect to derive a collisional life- 
time that is reasonably accurate with regards the evolution 
of the infra-red emission and of the planetesimal belt mass, 
but note that the collisional lifetime of objects of a specific 
size will not be quantitatively correct. In section [STT] we in- 
vestigate the effects of including a more realistic dispersal 
threshold that is dependent on size. 

Although the planetesimals in our model are not con- 
fined to a uniform ring for the entirety of the Nice model, we 
can still use this model to estimate the collisional lifetime 
by making the assumption that Rm is the radius contain- 
ing half of the mass of the disc, dr is the annuli containing 
98% of the mass and using the mean eccentricities and mean 
inclinations from the Nice model data. 

The mean radius of the belt is approximately constant 
at 26 AU during the pre-LHB phase and then rises dur- 
ing the LHB due to the planetesimals being scattered and 
reaches 79 AU. Similarly, the width of the belt is approx- 
imately 17 AU during the pre-LHB phase and rises to a 
maximum of 640 AU after the LHB. The radius and width 
of the belt at the end of the simulation are clearly much 
larger than the present day classical Kuiper belt. This is be- 
cause most of the objects left at the end of the simulation 
are scattered disc objects with a range of (typically high) 
eccentricities. Thus our assumption that the planetesimals 
are confined to a uniform ring clearly does not hold after the 
LHB, and moreover a range of collisional lifetimes would be 
expected depending on the objects' orbits. Nevertheless, we 
note that the collision time-scale derived above is within 
25% of that expected for an eccentric ring of planetesimals 
all with semi-major axes at 102 AU and eccentricities of 
0.56, which are the mean values from the simulation at the 
end of the post-LHB phase ([Wyatt et al]|2009[). Thus we 



expect the post-LHB collisional lifetimes presented here to 
be representative of an average member of the scattered disc 
in this phase, but use this simply to note that the collisional 
lifetime rapidly becomes larger than the age of the Solar 
System so that there is no further collisional mass loss (in 



agreement with Levison et al. 2008b I. A consideration of col 



lision rates in populations with a range of eccentricities and 
semi-major axes (see Wyatt et al.||20()9 \ would be required 
for a more detailed understanding of collision lifetimes in 
the post-LHB population. 

The mean eccentricities and inclinations give a mean 
relative velocity of around 360 m s~^ for the time before the 
LHB. At the time of the LHB, the large number of planetes- 
imals being scattered leads to a rapid increase in the mean 
relative velocity, which rises to a maximum of 2700 m s~^ 
at 900 Myr. After this time, the mean relative velocity grad- 
ually decreases to a value of 2500 m s~^ at the end of the 
simulation as the highly eccentric and inclined particles are 
more likely to be scattered out of the system. 
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Figure 7. Collision time-scale of the largest (Pluto-sized) objects 
as a function of time. Catastrophic collisions are only important 
before and during the LHB . The actual time-scales shown here are 
unrealistically small due to the assumption of a size-independent 
dispersal threshold. 



2.5.1 Collisional lifetime of the largest objects 



a bove, Qp should be m uch greater than 
Benz fc Asphaug||l999 1 for the largest ob- 



As discussed 

200 J kg-^ (e. 

jects (of size 2000 km), however, we can still use this model 
to estimate qualitatively the evolution of mass due to col- 
lisions. Figure [7| shows how the collision time-scale of the 
largest objects changes as the system evolves. The time-scale 
has been extrapolated assuming that the total mass is the 
only parameter in equation [9] that changes with time after 
the end of the Nice model simulation as described in section 
|2.2[ During the pre-LHB period, we find that the collision 
time-scale varies between 100-300 Myr. This implies that 
catastrophic collisions might have played a significant role 
in mass loss before the LHB. Since we need to end up with 
24 M0 of KBOs at the beginning of the LHB (see section 
2.11 we can approximate how massive the initial disc must 
have been to account for collisional mass loss. If we assume 
that before the LHB all the mass was lost through collisions 
then the total mass evolves as (e.g. Wyatt et al.|200fa k 

Mtot = M,„t/(1 4- t/t^D,)). (15) 

This means that we require an initial mass of ~150 to 
account for the mass lost due to collisions. By taking into ac- 
count the mass lost through dynamical processes (~10 M^) 
this gives us a rough estimate of 160 as the initial mass 
of the Kuiper belt, much greater than the 35 Me used in 



Gomes et al. (20051. Although this is a very rough approx- 



imation of the collisional evolution and these results are 
not quantitatively correct due to the assumption of a size- 
independent dispersal threshold (which will be considered 
in more detail in section [3| it does show us that collisions 
were important before the LHB and so would have affected 
the evolution of mass in the system. 

After the instability, the collision lifetime increases to 
beyond the lifetime of the Solar System, with a collision life- 
time of 130 Gyr by the end of the simulation and 5000 Gyr 
by the present day. This shows that collisions of the largest 
bodies become so infrequent that they can be neglected. 
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Figure 8. Collision, P-R drag and SW drag time-scales of the 
smallest particles (2.2 iiva) as a function of time. These time- 
scales assume that the particles are confined to a belt with a 
mean radius that increases due to the scattering caused by the 
LHB as described in the text. The P-R drag and SW drag time- 
scales become important after the LHB. 



2.5.2 Lifetimes of the smallest particles 



In section [T3l we assumed that the cut-off at the small end 
of the size distribution was defined by radiation pressure - 
particles smaller than Dbi will be blown out of the system by 
radiation pressure. However, particles larger than this may 
be affected by Poynting-Robertson (P-R) drag or solar wind 
(SW) drag on shorter time-scales than those for removal 
by collisions. To assess this we first calculate the collision 
time-scales using equation [9] for particles of size Dbi. This 
gives a time-scale that evolves similar to the time-scale for 
the largest objects but up to 5 orders of magnitude shorter 
(figure [8|. 

P-R drag is the tangential component of the radiation 
force which causes a decrease in both the semi-major axis 
and the eccentricity of a particle. The time-scale for a par- 
ticle of size Dhi to spiral into the Sun from a distance Rm 
under the influence of P-R drag is given by ( [Wyatt et al.| 
19991): 



tp_R = 400 



(16) 



where /3 — 0.5 for the smallest grains. 

A similar effect is caused by the tangential component 
of the solar wind known as the corpuscular stellar wind 
drag (hereafter SW drag). At the present time the SW drag 
force is equivalent to only 20-43% of the P-R drag force 
( Gustafson|1994 1, however, the higher mass loss rate of the 
young Sun means that SW drag would have been more effec- 



tive at removing dust than P-R drag at early times ( Minato 
et al.|2006[ ). The ratio of th e SW drag time-scale to the P-R 
drag time-scale is given by ( [Plavchan et al.|2005| ): 



tsw 

tp-R 



= 3.4- 



Mq 



(17) 



where Qp-r/Qsw is the ratio of the coupling coefficients, 
Mq is the present day mass loss of the Sun and M* is the 
mass loss of the Sun at different epochs. 



The History of the Solar System 's Debris Disc 9 



For this paper we have assumed that Qp-r/Qsw = 
1, Mq = 2 X 10~^* M0yr~^ and the luminosity remains 
constant at 1 Lq. Although the luminosity of the Sun has 
changed with time and is likely to have been only 0.7 L0 



when it first became a main-sequence star (e.g. Jorgensen 



1991 and references therein), the uncertainties in the other 
factors in this equation are much greater than those due to 
this change in luminosity. For the change in stellar mass loss 
rate with time we have taken M*(t) = Af0(f/4.5 Gyr)"^-^^ 
fo r t > 700 Myr from the analysis of stellar mass loss rates 



Wood et al. 



(20051 and M*(t) =80 Mq for earlier times. 



Figure [8] compares the collision time-scale with the P- 
R drag time-scale and the SW drag time-scale for parti- 
cles of size Dhi- CoUisional processes would have dominated 
the removal of dust before the LHB due to the high dust 
mass present and the compactness of the belt. Drag forces 
would have been insignificant as they are in all other ob- 
served debris discs ( |Wyatt|2 005). However, as the distribu- 
tion of mass becomes increasingly spread out during and af- 
ter the LHB, the collision time-scales of the particles rapidly 
increases such that they become more susceptible to drag 
forces. Due to the high mass loss rate of the early Sun, SW 
drag is more effective at removing dust than P-R drag until 
~2.7 Gyr. The increased rate of dust removal due to drag 
forces throughout the post-LHB phase reduces the amount 
of small dust below that expected in the collisional cas- 
cade equation |10[ This means that the collisional lifetime 
of the smallest particles is in fact underestimated in figure [S] 
(which assumes the collisional cascade size distribution ex- 
tends down to the blow-out limit). Since it is the smallest 
dust which contributes most to the emission, this increased 
rate of dust removal also reduces the emission and means 
that we have overestimated the fractional luminosity (shown 
in figure [4|) and the observable properties dependent on this 
(figures [5|and[6]l at late times. However, since it is predicted 
that the Kuiper belt would not be observable at this time 
(figures [5] and [6| , this does not affect the comparison with 
observed debris discs. 



We note that the time-scales used here are all for the 
average particles in the model. In the pre-LHB phase, some 
particles are occasionally scattered in from the narrow belt 
making them more susceptible to drag forces. During and 
after the LHB, the range of orbital elements of the particles 
is greatly increased. Since both of the drag forces are pro- 
portional to R^, the position of a particle will greatly affect 
whether it is destroyed through collisions or spirals into the 
Sun due to drag forces. 

We also note that the mass loss rate of the Sun at early 



times is not very well constrained. Some authors (e.g. Sack- 
|mann fc Boothroyd|2003l ) have suggested that the mass loss 
rate of the Sun in the early Solar System may have been as 
much as 1000 times greater than the current value. If the 
mass loss rate was this high then the SW drag time-scale 
would have been shorter than the collisional time-scale for 
the smallest particles reducing the small size end of the size 
distribution, thus reducing the luminosity of the disc in the 
pre-LHB phase below that presented here. 
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Figure 9. Evolution of the size distribution for an initial transi- 
tion diameter of 70 km and initial mass of 40 Mgj. The size dis- 
tribution continues down to D]-,i but this has not been shown for 
clarity. The thick dashed lines represent the break diameter and 
the transition diameter. There is some increase in Dt pre-LHB 
when mass is being lost through collisions but the size distribu- 
tion then remains fixed from the onset of the LHB and all the 
mass lost after this time is due to dynamics. 



3 REALISTIC SIZE DISTRIBUTION AND 
GRAIN PROPERTIES 

In section [2] we described a basic method for investigating 
the history of the Solar System's debris disc. That modelling 
method can readily be applied to the outcome of any numer- 
ical simulation to consider its observable properties. In this 
section we will confront two of the main assumptions of the 
model. So far we have been assuming that the planetesimals 
and dust are governed by a single phase size distribution 
and that they are black bodies (with a slight correction at 
sub-mm wavelengths). 



3.1 Three-phase size distribution 



In section [2.3| we assumed that the particles are in a colli- 
sional equilibrium from the smallest to the largest particle. 



Using their collisional evolution code, Lohne et al. ( 2008 1 



show that a disc that starts with a single power law size 
distribution will quickly develop into a system with a three- 
phase power law due to differences in the collisional time- 
scales of different sized particles. As the system evolves (fig- 
ure |9|, the particles begin to reach collisional equilibrium 
starting with the smallest particles due to their shorter col- 



lisional lifetime (see section 2.5 1. The transition diameter, 
Dt, defines the diameter at which the collisional lifetime of 
the particles is equal to the age of the system. Particles be- 
low this size will reach collisional equilibrium. The slope of 
the power law for particles smaller than Dt will then de- 
pend on whether the particles are in the strength or gravity 
regimes. The slope of the power law for particles larger than 
Dt is given by the primordial slope, qp. 

In section [2^ we assumed that the dispersal threshold 
of a particle is independent of its size. This dispersal thresh- 
old defines the minimum energy required to catastrophically 
destroy a planetesimal and disperse the fragments such that 
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they do not recombine under their own gravity and is, in 
fact, dependent on the size of the particle. The dispersal 
threshold decreases with size for the smallest size particles 
for which little energy is required to disperse the fragments 
after the collision (the strength regime, see e.g. |Farinella 



eFal] [19821 |Housen fc Holsapple| [19901 |Benz fc Asphaug 
1999). As object size increases so does the gravitational 



strength (assuming constant density) to the extent where 
the energy required to disperse the fragments of a collision is 
greater than the energy required to catastrophically destroy 
the planetesimal. For these objects, the dispersal threshold 



increases as size increases (the gravity regime, see e.g. Pe- 



tit fc FarineIIa|[T993| |Campo Bagatin et al.[[T994| |Benz fc 



Asphaug|1999 1. These two regimes can be described by the 
sum of two power laws (e.g. Krivov et al.||2005 l: 



>d(^) = A 



D 

27n 



D 
2 km 



(18) 



where A^, A^, bs and 6g are free parameters. This equation 
can then be used to find the transition between the strength 
and gravity regimes (the diameter at which the two power 
law components contribute equally), which occurs at the 
breaking diameter, Db (in km): 

l/{3bg-3(),) 



As 



(19) 



Ag (2 X 10-3)36= 

To find Dt we assume that the disc starts as a single 
power law with a primordial distribution given by qp. Using 
equation [9] (and adjusting atot and fee for a three-phase 
distribution) we can find the diameter for which an object's 
coUisional time-scale is the same as the age of the system, 
tc(-Dt) = t. Objects smaller than Dt will be governed by 
either the strength regime (if their size is also below Dh) 
with a slope of qs or the gravity regime with a slope of q^. 

If we know the transition diameter at a particular time 
then we can work out the total mass of the disc at that time: 



Mtoti(t) = 
B{t) = 



8.7 X W-"n,^^4t)pB{t) 

Wmax(O) 

iTt/MAj 



(20) 
(21) 



A(t) 



3'?g-3gp£)393~3ij, 



6 — 3gs 



6 -39b 



D oqg 



d: 



6 — 3(3p 



6 — 3gp 



(22) 



where rimax is the number of objects of size Dc and assuming 
that none of qa, q^ and qp are equal to 2. Given the param- 
eters of our model, we find that the largest objects can only 
be destroyed by objects much larger than themselves and so 
cannot be coUisionally destroyed in any of our simulations 
(see section p{.1.2 1 allowing us to ignore the t/tc{Dc) term. 

However, since Dt{t) is dependent on Mtoti(i) (through 
equations [9l and [To| , Dt is calculated at each time-step based 
on the total mass of the previous time-step. If we ignore the 
dynamical mass loss from the system, then the total mass 
at each time-step is given by: 



Mtotl(U) = Mtotl(tn-l) 



B{tn 



(23) 



If we assume that the coUisional evolution does not affect 
the rate of mass lost from dynamical evolution then the dy- 
namical evolution of the Nice model (as described in section 
2.2 I can be combined with this equation to give us a total 
mass, Mtoti, that evolves as: 



Mtotl(tn) = Mtotl(tn- 



B(t„) 



AMtot(tn) 



(24) 



^B(t„-l) Mtot{tn-l) , 
AMtot(tn) = Mtot(tn) - Aftot(tn-l). (25) 

By keeping our assumption that the cross-sectional area 
is proportional to the mass and that the proportionality con- 
stant is defined by the size distribution, the changing size 
distribution (figure [9| can then be used to find how the ra- 
tio of cross-sectional area to mass changes with time. This 
cross-sectional area can then be used to find the emitted flux 
as described for the basic model in section l273l 



3.1.1 Parameter choices 

The free parameters in this model are Qd, gp, -Dt(O) and 
Aftoti(O). Lohne et al. (20081 use a Qd similar to that 
of Benz & Asphaug (19991. For the coefficients they set 
As = Ag = 500 J kg- and the exponents are 36s ~ —0.3 and 
36g = 1.5. The transition between the strength and gravity 
regimes (the diameter at which the two power law compo- 
nents contribute equally) occurs at the breaking diameter, 

Db = 632 m. 

Leinhardt & Stewart (20091 show that, if comets have 
negligible strength, then their Qn function can be much 
lower than those of Lohne et al. ( 2008 1 and Benz fc As- 



phaug ( 1999 1, with coefficients as low as A, 
^ and exponents of 36a 



Ag = 28 J kg 
which gives a breaking diameter of Db = 324 m 



20 J kg-^ and 
0.4 and 36g = 1.3, 



Stewart & Leinhardt ( 2009 \ go on to show that using 



Qd is only valid when the target object is much bigger 
than the object impacting it. When the impacting object 
is roughly half the size of the target object or larger (assum- 
ing impact and target have equal density), Qq is no longer 
valid as it only takes into account the size of the target ob- 
ject rather than the size of both objects. In our simulations, 
this means that objects of size >130 km become harder to 
catastrophically destroy and objects >440 km become im- 
possible to destroy. However, as the objects most strongly 
affected by this are bigger than the final transition diam- 
eter we are interested in, this does not make a significant 
difference to our work. 

As the slope of the primordial distribution remains con- 
stant with time, this parameter can be found from present 
day observations. Recent surveys of the largest KBOs give 
a size distribution slope qp — 1.8 — 2.3 (see Petit et al 
and references therein). Although Bernstein et al. 



2008 



(2004h 



show that the classical belt and the excited belt have differ- 
ent size distributions, for this work we shall consider both 
populations to have the same size distribution. The most re- 
cent observations suggest a slope of gp 2.2 ( Fuentes et al. 
2009[|Fraser fc Kavelaars|2009| ) and so we will use this value 
in the rest of this paper, qs and Qp can be found from the 



formula ( [O'Brien fc Greenberg|2003 1 

_ 22/6 6 



(26) 
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which sets qs = 1.877 and Qg = 5/3. 

The initial conditions, -Dt(O) and Mtoti(O) 
([2008) 



200 



constrained. Lohne et al 



are less well 
start their model with a 
single size distribution and sec how it evolves in to a three- 
phase size distribution. As the smallest particles have such 
a short collisional lifetime, the size distribution is likely to 
have already undergone some evolution before the time at 
which the Nice model begins and so -Dt(O) 2> -Dbi. Numerical 
simulations of the accretion of KBOs suggests that the initial 
transition diameter should be no less than 1 m and probably 



more than 100 m (Kenyon & Luu 1999 1. For the asteroid belt 



it has been shown that the transition diameter is Ukely fixed 
during the accretion phase ( Bottke et al.|2005 l. If the same 
is true of the Kuiper belt then we would expect an initial 
transition diameter of ~ 100 km. Models of formation of 
Kuiper belt objects require that there must have been at 
least 10 in the primordial Kuiper belt for the largest 
KBOs to have formed (e.g. |Stern||1996a| ). 

However, there are stronger constraints on the values 
of Dt and Aftoti at the time of the LHB. The Nice model 
requires that there must be 24 of planetesimals in the 
disc for the LHB to occur and that the size distribution be- 
comes fixed at this time and so the transition diameter must 
be equal to the transition diameter in the current Kuiper 
belt. Observations constrain this diameter to between 50 and 
200 km (Ber nstein et al.|2004| [Fuentes et al.|2009[|Fraser fc 
IKavelaars 20091. Taking these constraints into account, we 



can run the simulation with various initial conditions and 
find out which give the appropriate results. 

The top plot of figure [lO] shows simulation runs using 
the Lohne et al. ( 2008 ) . The white region shows the runs 



which give a final transition diameter within the constraints 
of current observations and the 24 line shows the runs 
which leave us with enough mass in the planetesimal belt for 
the LHB to occur. From this we can see that if we start with 
a transition diameter >100 km then there is no collisional 
mass loss since objects of this size have a collision timescale 
longer than 879 Myr. 

The bottom plot of figure 10] sh ows simulation runs us- 
ing the Leinhardt & Stewart (20091 Q^. In these runs the 
24 line does not overlap with the white area showing 
that we cannot satisfy both of our constraints with this Q^. 



This means that the Leinhardt & Stewart ( 2009 1 is in- 



consistent with the Nice model interpretation of the LHB as 
not enough mass remains in the Kuiper belt for the required 
length of time. This may be a result of our assumption that 
the dynamical evolution is not affected by the collisional 
evolution; it could mean that the objects are stronger than 
one would expect; or it might actually mean that the Nice 
model cannot be used to explain the LHB. To remain consis- 
tent with the Nice model interpretation of the LHB, we shall 
use the stronger QJj throughout the rest of this paper. We 
will also set qp = 2.2, _Dt(0) = 70 km and Mtoti(O) = 40 M® 
to provide an illustrative case where the mass at the time 
of the LHB is 24 M® and the final transition diameter is 
87 km. 



3.1.2 Effect of a three-phase size distribution 

In section [2] we used a size distribution with a slope defined 
by qd = 11/6 to describe all of the objects. This results in 
a distribution where most of the mass is concentrated in 
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Figure 10. Dependence of mass at the time of the LHB on the 
initial conditions. Contours show the mass in Earth masses and 
the shaded regions represent runs that give a final transition diam- 
eter smaller or larger than the current constraints on the present 
day transition diameter. The runs shown in the top plot use the 
[Lohne et al.| l|2008[l and qp = 2.2 and those in the bottom 
plot use the [Leinhardt fc Stewart] | |2009| and qp = 2.2. 



the largest objects and so the collisional mass loss of the 
system was dependent on the collisional time-scale of the 
largest objects. We now find that the largest objects of size 
2000 km now require a dispersal energy of Qy){Dc) — 1.6 x 
lO'^ J kg~^ rather than the 200 J kg~^ used in section [2] 
This now means the largest objects can only be destroyed 
by objects much larger than themselves and hence will never 
be destroyed in our simulations. Therefore, the largest KBOs 
in the Solar System are likely to be primordial as also found 



by the work of Farinella & Davis (19961. However, setting 
the slope of the primordial size distribution to 2.2 means 
that most of the mass is concentrated in objects with sizes 
7? « Dt. As -Dt increases, mass is still lost through collisions 
as the primordial planetesimals reach collisional equilibrium. 

Figure [TT] shows how the SED changes when we use a 
3-phase size distribution. By comparing the lines for black- 
body 1-phase and black-body 3-phase we can see that the 
3-phase model increases the amount of flux emitted at all 
wavelengths, with a 4-fold increase in the peak of the SED. 
Figure [12] shows how the evolution of /, F^/Fv* and Mdust 
changes when we use a 3-phase size distribution. From this 
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Figure 11. SED pre-LHB for black-body grains with a single 
slope size distribution and also black-body grains, amorphous sil- 
icate grains and comet-like grains with the Lohne model of an 
evolving size distribution. The realistic grain models also include 
a three-phase size distribution and all of the three-phase models 
have an initial mass of 40 M0 and an initial transition diameter 
of 70 km. 



plot we can see that the increase in flux is present at all 
times. This increase in flux is due to the increase in the 
cr/M ratio. At early times it is five times greater than that 
given in equation [2] and decreases to four times greater as 
Dt increases. This slight decrease can be seen in the dotted 
lines of figure [12] where the slope of the line becomes slightly 
steeper just prior to the LHB. 

This result implies that before the LHB, the Solar Sys- 
tem would have been amongst the brightest debris disc sys- 
tems around Sun-like stars in 24 /xm emission and 70 /im 
emission. Although plausible, this could be a result of our 
choice of parameters. For instance, in equation [18] we set 
feg = 0.5 giving us a power law slope in the gravity regime 
of = 5/3, whereas recent observations suggest a value of 
gg « 4/3 might be more realistic ( Fuentes et al.|2009 Fraser 
fc Kavelaars|2009[ ). 

By changing from a single slope power law to a three- 
phase power law and including a dispersal threshold that is 
dependent on size means that the Solar System's debris disc 
would have been significantly brighter, but the rest of the 
conclusions of section |2] still hold. The 24 /^m and 70 /im 
excesses still evolve in the same manner with a peak in the 
24 jim. excess still being seen at the time of the LHB. Both 
models also show that mass lost from collisions is likely to be 
important in the evolution of the Solar System's debris disc, 
although in the basic model the mass loss was from collisions 
of the largest objects whereas now we have shown that the 
mass loss is due to the evolution of the size distribution and 
the amount of mass lost is greatly dependent on the initial 
size distribution used. 



3.2 SED model 

In reality, particles do not emit or absorb efficiently at all 
wavelengths. Their emission and absorption efficiencies are 
dependent on their composition and porosity. | Wyatt fc Dent| 
( 2002 1 created a model that shows how the SED changes 



compositional model of Li & Greenberg (19971. They used 



this model to find the most likely composition of the Fo- 
malhaut debris disc. Here we use this model to see how it 
changes the results of section |2] 

By relaxing the black-body assumption we change equa- 
tions [3] and |4] so that: 



= 2.35 X 10-"d-'^^a(i?)Qabs(A,D)r7(D) 



T{D,R) 



(27) 
(28) 



where Tbb is the black-body temperature given by equa- 
tion [4] and Qabs (the absorption efficiency) is found using 
Mie theory, Rayleigh-Gans theory or geometric optics in the 
appropriate limits using optical properties from the compo- 
sitional model. 

To calculate the composition of the particles, the core- 



mantle model of Li & Greenberg (19971 is used. This as- 



sumes the particles to be formed of a silicate core surrounded 
by an organic refractory mantle where the silicate core makes 
up 1/3 of the total volume. The silicate material may be ei- 
ther amorphous or crystalline. The particles will have a given 
porosity (p) defining how much of the particles' volume is 
empty space, and a given fraction of water (gH2o) defining 
how much of the empty space is filled by ice. 

Although the composition of some of the largest KBOs 
has been inferred from spectroscopy (e.g. Barucci et al. 



2008 1 , it is the smaller grains (which have not been ob- 



served) that have a larger effect on the SED. In this paper 
we will look at two extremes for the particle composition. 
The first grain composition replaces the black-body grains 
by amorphous silicate grains with zero porosity and no ice. 
For this composition p = 2370 kg m~^ and Dhi = 1.47 /xm. 
For the second grain composition we will assume that the 
grains have a similar composition to the "comet-like" grains 
of Augereau et al. (19991. As such, this composition uses 
crystalline grains with p = 0.93 and gH2 = 0.38. For this 
composition p = 590 kg m~'^ and Dbi = 0.85 pm. 



3. 2. 1 Effect of using a realistic grain model 

Figure [11] shows how the SED for the realistic grain models 
compares to the black-body models from sections [2] and [3d" 
(both realistic grain models also use the evolving size distri- 
bution described in section 



3.1 1. From this plot we can see 



that by introducing realistic grains to the model, the peak 
of the SED clearly moves to a lower wavelength. This oc- 
curs because the majority of the emission is from particles 
in the size regime where they absorb radiation more effi- 
ciently than they emit it, which causes their temperature to 
increase more than if they were black-bodies (see equation 
28 1. The SED also shows features in the emission due to the 
elements present in the grain, although these are much more 
prominent in the "comet-like" model due to the presence of 
ice and the crystalline nature of the grains. 

The fact that the peak of the wavelength emission has 
moved to a lower wavelength means that the 24 pm flux is 
much higher compared to the black-body grains (see figure 
12 L So high that the Solar System would have been amongst 
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Figure 12. Same as figures [4)61 but also including the black-body three-phase size distribution model and the realistic grain models. 
The realistic grain models also include a three-phase size distribution and all of the three-phase models have an initial mass of 40 IVI^ 
and an initial transition diameter of 70 km. 



the brightest discs at 24 fim before the LHB. This increase 
in 24 jj,m flux at young ages means that we no longer see 
the jump in ^24/^24* that we saw for the black-body grains 
at the time of the LHB. However, this jump is seen at lower 
wavelengths, such as those below 10 ^m, as can be seen in 
figure [13] which shows the SED for the "comet-like" grains 
before, during and after the LHB. The 70 fim flux is reduced 
when using realistic grain properties, although this is more 
pronounced for the amorphous silicate grains. This decrease 
is because the SED peak is at a lower wavelength and the 
total flux has decreased due to the decrease in emitting ef- 
flciency of the particles. 

Using realistic grain models changes the density from 
that assumed previously, especially when porous grains are 
used. The lower density of the "comet-like" grains has the ef- 
fect of increasing the a/M ratio which reduces the coUisional 
lifetime of all particles. This reduction in coUisional lifetime 
means that a greater flnal transition diameter is reached. 
For the case of the "comet-like" grains, a flnal transition di- 
ameter of 117 km is reached. However, it should be noted 



that we have assumed objects of all sizes to have the same 
density and porosity, whereas the large KBOs have lower 
porosities and so higher densities - Varuna, for example, is 



estimated to have a porosity of 0.05-0.3 ( Jewitt & Sheppard 



2002 1 - and therefore there would not be as much coUisional 



evolution as we have suggested here. If amorphous silicate 
properties are used, the density of the planetesimals is much 
higher and so they do not collisionally evolve and the tran- 
sition diameter does not change. 

The high water-ice content of the "comet-like" model 
means that these particles are likely to undergo ice sublima- 
tion close to the Sun, however, the precise effects of subli- 
mation on the size distribution are unknown. Since particles 
would only come close enough to the Sun for a brief period 
during the LHB (see figure [T| , we would only expect ice 
sublimation to affect our results during the brief mid-LHB 
phase and we intend to explore this effect in a future work. 

Introducing realistic grains to our model has the effect 
of moving the SED to shorter wavelengths, thus increasing 
the mid-IR flux emitted. The same effect could be approx- 
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Figure 13. Same as figure [3] but for the comet grain model. 

imated with black-body grains if they were assumed to be 
much closer to the star and therefore have a greater temper- 
ature. 

Of the four models presented in this paper, the "comet- 
like" is the most realistic. This predicts that the current 
fractional luminosity is 3 x 10~^ and the current submil- 
limetre dust mass is 7 x 10~^ Mm. Both of which are within 



the observational limits (see sections 2.3 and 2.4 1 



4 CONCLUSIONS 

In this paper we present a new look at the history of the So- 
lar System. Starting with the Nice model for the evolution of 
the Solar System, we demonstrate how the evolving spatial 
distribution of planetesimals in the system changes the ther- 
mal emission of dust produced in planetesimal collisions. We 
started with a simple model that used a single-phase power 
law to convert between mass of planetesimals and surface 
area and assumed black-body emission (see section |2|. We 
find that this model predicts that the primordial Kuiper 
belt would have been detectable at 24 and 70 /xm before the 
LHB. During the LHB more hot dust would have been pro- 
duced as many of the KBOs are scattered inwards towards 
the Sun causing a peak in the 24 /im emission. Within a few 
hundred million years of the LHB, the dynamical depletion 
of the Kuiper belt renders it undetectable at 24 and 70 /xm. 



Statistics from surveys of Sun-like stars (e.g. Trilling 
et al. 2008 1 show that the number of stars with an observable 



excess at 24 ^ra decreases with stellar age and the number 
of stars with an observable excess at 70 ^m remains approx- 
imately constant with stellar age. An LHB-like event causes 
a drop in excess ratio of approximately 4 orders of magni- 



tude at both these wavelengths (see figure 12 1 showing that 
such major clearing events must be rare and that most de- 
bris discs that are detectable just after 10 Myr lose mass 
through collisions rather than through dynamical instabil- 
ities. This allows us to set an upper limit of 12% on the 
fraction of Sun-like stars that go through an LHB. 

Figure [5] shows that the period of increased 24 /xm emis- 
sion only lasts for ~15 Myr. If we assume that the average 
age of Sun-like stars that we observe is 5 Gyr then there is 
at most a 0.04% chance of observing a star going through 



a LHB. However, the bombardment from the asteroids lasts 
about 5 times longer ( Gomes et al.|[2005 l and so could in- 
crease this possibility to 0.2%. Therefore, of the 26 Sun-like, 
field stars older than 10 Myr found to have a 24 ^m ex- 



cess out of a compiled sample of 413 (see table 5 in Gaspar 
et al. 20091, approximately 1 could be an observation of a 



current LHB event. Certain systems such as rj Corvi which 
have both hot and cold dust but too much hot dust to be 
explained by coUisional evolution alone ( Wyatt et al.|2007a l 
may still be explained by an LHB-like event. 

Although in this model we have just considered the evo- 
lution of the Solar System as described by the Nice model, 
this paper gives enough detail for this simple model to eas- 
ily be applied to the output of any numerical simulation of 
planetesimals to give an indication of its observable proper- 
ties. 

In section |3] we removed the major assumptions of the 
initial model by including a three-phase power law size dis- 
tribution that depends on coUisional history. Changing the 
size distribution has the effect of greatly increasing the flux 
emitted. Having a size distribution that changes during the 
simulation also means that mass is also lost due to collisions. 
For instance, if we start with an initial transition diameter of 
70 km, an initial mass of 40 M® is required to leave us with 
the 24 Mq at the time of the LHB that is necessary for the 
LHB to take place. However, this assumes that the KBOs 
are as strong the Benz &: Asphaug|((1999 ) case. If they are as 
weak as Leinhardt & Stewart ( 2009 1 suggest (and our sim- 
plistic combination of coUisional and dynamical mass loss 
is correct) then this shows that the Nice model cannot be 
used to explain the LHB since either too much mass is lost 
through coUisional grinding or the final transition diameter 
is unrealistic. This does not rule out the Nice model in its 
entirety only the need for there to be a delay before the 2:1 
mean motion resonance crossing of Jupiter and Saturn. 

We also find that changing the grain properties to re- 
semble more realistic grains has the effect that the spike in 
24 /im emission is no longer seen as the peak wavelength is 
shorter and so the 24 yLxn emission is initially much higher. 
However, peaks in emission at lower wavelengths such as 
10 /im would still be possible indicators of an LHB-like, 
transient event. The wavelengths the spike appears at also 
depends on the heliocentric distance of the disc. If the disc 
starts further out in the system then it could still give a 
spike in the 24 ^m emission. 

One major caveat of this work is that we have only 
concentrated on the emission from the Kuiper Belt. Since 
this is far from the Sun most of the emission will be from cold 
dust and thus we are underestimating the warm emission 
and, therefore, the 24 /im fiux. In future work we intend to 
include the asteroid belt in this model and investigate the 
effect of sublimation from comets. 
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